
setwd("..")
setwd("..")
setwd("..")
setwd("..")

Dawson_Path <- "Project_DisadvantagedCommunities/Code/CVRP_EFMP_EXP_repo/Dawson/"

library(tidyverse)
library(lubridate)
library(stringr)

# Mapping
library("raster") # Shape File Manipulation
library("rgeos") # Geospatial
library("rgdal") # Geospatial
library("ggmap") # Google Maps Layers

archive <- Sys.Date() %>%
  ymd() %>%
  str_replace_all(., "-", "")

theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_blank(),
                         plot.background = element_rect(fill=NA),
                         panel.border = element_blank(),
                         axis.line=element_blank(),
                         axis.text.x=element_blank(),
                         axis.text.y=element_blank(),
                         axis.ticks=element_blank(),
                         axis.title.x=element_blank(),
                         axis.title.y=element_blank(),
                         plot.title = element_text(size=22)))

## Opening Data
# DAC Status
test_data <- readstata13::read.dta13("Project_DisadvantagedCommunities/Code/CVRP_EFMP_EXP_repo/Dawson/test.dta") %>%
  as_tibble()
# CA Outline
states <- readOGR(dsn=paste0(Dawson_Path, "data/states_21basic"), layer="states")
CA <- states[states$STATE_ABBR %in% "CA", ] # Isolate CA
CA_df <- fortify(CA); rm("states")

CAcities <- map_data("state", region="california")
CAcities <- data.frame(State=rep("california",9), 
                       City=c("Los Angeles", "San Diego", "San Francisco", 
                              "Sacramento", "Modesto", "Bakersfield",
                              "San Jose", "Fresno", "San Luis Obispo"))
CAcities <- cbind(geocode(as.character(CAcities$City)), CAcities)

air_basins <- readOGR(dsn=paste0(Dawson_Path, "data/ca_air_basins"), layer="CaAirBasin") %>%
  spTransform(., crs(CA))
air_basins@data$id <- rownames(air_basins@data)
air_basins_df <- fortify(air_basins) %>%
  left_join(., air_basins@data)
air_basins_pilot_df <- air_basins_df %>%
  filter(NAME %in% c("South Coast", "San Joaquin Valley"))

air_basin_names <- aggregate(cbind(long, lat) ~ NAME, data=air_basins_pilot_df, 
                             FUN=function(x)mean(range(x)))
air_basin_names$long <- air_basin_names$long + c(0.33, 1.1)
air_basin_names$lat <- air_basin_names$lat + c(-0.1, 0)

ZIP_poly <- readOGR(dsn=paste0(Dawson_Path, "data/cb_2016_us_zcta510_500k"), layer="cb_2016_us_zcta510_500k")

## Adding DAC status
DAC_Status <- test_data %>% distinct(zip, dac)
ZIP_poly@data$id <- rownames(ZIP_poly@data)
ZIP_poly <- ZIP_poly[CA,]
ZIP_poly_df <- fortify(ZIP_poly) %>%
  left_join(., ZIP_poly@data)
ZIP_poly_df <- ZIP_poly_df %>%
  left_join(DAC_Status, by = c("GEOID10" = "zip"))

########
## Mapping
########
ggplot() +
  geom_polygon(data=CA_df, aes(long,lat, group=group), color="black", size=0.4, fill = NA) + 
  geom_point(data = CAcities, aes(x = lon, y = lat), size = 0.25) + 
  geom_text(data = CAcities, aes(x = lon, y = lat, label = City), 
            hjust = 0, vjust = 0, nudge_x = 0.1, size = 3) +
  theme_opts
ggsave(paste0(Dawson_Path, "maps/map_CA_outline_",archive, ".png"))

# CA with Air Basins
ggplot() +
  geom_polygon(data=CA_df, aes(long,lat, group=group), color="black", size=0.4, fill = NA) + 
  geom_polygon(data=air_basins_df, aes(long,lat, group=group), color="black", size=0.3, fill = NA) +
  geom_point(data = CAcities, aes(x = lon, y = lat), size = 0.25) + 
  geom_text(data = CAcities, aes(x = lon, y = lat, label = City), 
            hjust = 0, vjust = 0, nudge_x = 0.1, size = 3) +
  theme_opts
ggsave(paste0(Dawson_Path, "maps/map_AQMD_Outline_",archive, ".png"))

# CA with Air Basins and Pilots Shaded
ggplot() +
  geom_polygon(data=CA_df, aes(long,lat, group=group), color="black", size=0.4, fill = NA) + 
  geom_polygon(data=air_basins_df, aes(long,lat, group=group), color="black", size=0.3, fill = NA) +
  geom_polygon(data = air_basins_pilot_df, aes(long, lat, group=group), color="black", size=0.3, fill = "grey") +
  geom_text(data = air_basin_names, aes(long, lat, label = NAME), fontface = "bold") +
  geom_point(data = CAcities, aes(x = lon, y = lat), size = 0.25) + 
  geom_text(data = CAcities, aes(x = lon, y = lat, label = City), 
            hjust = 0, vjust = 0, nudge_x = 0.1, size = 3) +
  theme_opts
ggsave(paste0(Dawson_Path, "maps/map_AQMD_Pilot_Fill_",archive, ".png"))

## Zoom in
ggplot() +
  geom_polygon(data=CA_df, aes(long,lat, group=group), color="black", size=0.4, fill = NA) + 
  geom_polygon(data=air_basins_df, aes(long,lat, group=group), color="black", size=0.3, fill = NA) +
  geom_polygon(data = air_basins_pilot_df, aes(long, lat, group=group), color="black", size=0.3, fill = "grey") +
  geom_text(data = air_basin_names, aes(long, lat, label = NAME), fontface = "bold") +
  geom_point(data = CAcities, aes(x = lon, y = lat), size = 0.25) + 
  geom_text(data = CAcities, aes(x = lon, y = lat, label = City), 
            hjust = 0, vjust = 0, nudge_x = 0.1, size = 3) +
  coord_map(ylim =  c(min(CA_df$lat), 38.5)) +
  theme_opts
ggsave(paste0(Dawson_Path, "maps/map_AQMD_Zoom_",archive, ".png"))

## Add ZIPs
ggplot() +
  geom_polygon(data=CA_df, aes(long,lat, group=group), color="black", size=0.4, fill = NA) + 
  geom_polygon(data=air_basins_df, aes(long,lat, group=group), color="black", size=0.3, fill = NA) +
  geom_polygon(data = air_basins_pilot_df, aes(long, lat, group=group), color="black", size=0.3, fill = "grey") +
  geom_polygon(data=ZIP_poly_df, aes(long,lat, group=group), color="black", size=0.2, fill = NA) +
  theme_opts + coord_map(ylim =  c(min(CA_df$lat), 38.5))
ggsave(paste0(Dawson_Path, "maps/map_Zoom_ZIPs_",archive, ".png"))

## Add ZIPs and fill in DACs
ggplot() +
  geom_polygon(data=CA_df, aes(long,lat, group=group), color="black", size=0.4, fill = NA) + 
  geom_polygon(data=air_basins_df, aes(long,lat, group=group), color="black", size=0.3, fill = NA) +
  geom_polygon(data = air_basins_pilot_df, aes(long, lat, group=group), color="black", size=0.6, fill = "grey") +
  geom_polygon(data=ZIP_poly_df, aes(long,lat, group=group), color="black", size=0.2, fill = NA) +
  geom_polygon(data = ZIP_poly_df %>% filter(dac == 1), aes(long, lat, group=group), color="black", size=0.3, fill = "red", alpha = 0.3) +
  theme_opts + coord_map(ylim =  c(min(CA_df$lat), 38.5))
ggsave(paste0(Dawson_Path, "maps/map_Zoom_ZIPs_DACs_", archive, ".png"))

## Zooming into the DACs
hold_Zoom <- ggplot() +
  geom_polygon(data=CA_df, aes(long,lat, group=group), color="black", size=0.4, fill = NA) + 
  geom_polygon(data=air_basins_df, aes(long,lat, group=group, fill = NA), color="black", size=0.2, fill = NA) +
  geom_polygon(data = air_basins_pilot_df, aes(long, lat, group=group, fill = "grey"), color="black", size=0.7, fill = "grey") +
  geom_polygon(data=ZIP_poly_df, aes(long,lat, group=group), color="black", size=0.2, fill = NA) +
  geom_polygon(data = ZIP_poly_df %>% filter(dac == 1), aes(long, lat, group=group, fill = "red"), color="black", size=0.2, fill = "red", alpha = 0.3) +
  scale_fill_manual(values = c(NA, 'grey', 'red'),
                    labels = c('Non-AQMD/Non-DAC', 'AQMD', 'DAC')) + 
  theme_opts + theme()
## Zoom SF
hold_Zoom + coord_map(xlim = c(-123, -121), ylim =  c(37, 38.5))
ggsave(paste0(Dawson_Path, "maps/map_Zoom_ZIPs_DACs_SF_",archive, ".png"))
## Zoom LA
hold_Zoom + coord_map(xlim = c(-119, -117.5), ylim =  c(33.5, 34.5))
ggsave(paste0(Dawson_Path, "maps/map_Zoom_ZIPs_DACs_LA_",archive, ".png"))
## Zoom San Jose
hold_Zoom + coord_map(xlim = c(-122.5, -120.5), ylim =  c(36.75, 38))
ggsave(paste0(Dawson_Path, "maps/map_Zoom_ZIPs_DACs_SanJose_",archive, ".png"))
## Zoom Sacremento
hold_Zoom + coord_map(xlim = c(-122.5, -120.5), ylim =  c(37.5, 39))
ggsave(paste0(Dawson_Path, "maps/map_Zoom_ZIPs_DACs_Sac_",archive, ".png"))









